!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!>\brief
!!
!! Volcano eruption.
!!
!! \n
!!
!<----------------------------------------------------------------------
module mod_volcano_test

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_fu_manager, only: &
   mod_fu_manager_initialized, &
   new_file_unit

 use mod_physical_constants, only: &
   mod_physical_constants_initialized, &
   t_phc

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_volcano_test_constructor, &
   mod_volcano_test_destructor,  &
   mod_volcano_test_initialized, &
   test_name,        &
   test_description, &
   phc,              &
   ntrcs,            &
   ref_prof, t_ref,  &
   coeff_dir,        &
   coeff_neu,        &
   !coeff_sponge,     &
   coeff_visc,       &
   coeff_init,       &
   coeff_outref

 private

!-----------------------------------------------------------------------

 ! public interface
 character(len=*), parameter   :: test_name = "volcano"
 character(len=100), protected :: test_description(3)
 type(t_phc), save, protected  :: phc
 integer, parameter            :: ntrcs = 0
 character(len=*), parameter   :: ref_prof = "isothermal"
 real(wp), target              :: t_ref

 logical, protected :: mod_volcano_test_initialized = .false.
 ! private members
 real(wp) :: &
   nu,     & ! kinematic viscosity
   p_vent, & ! vent pressure
   t_vent, & ! vent temperature
   v_vent    ! vent velocity
 character(len=*), parameter :: &
   test_input_file_name = 'volcano.in'
 character(len=*), parameter :: &
   this_mod_name = 'mod_volcano_test'

 ! Input namelist
 namelist /input/ &
   nu, p_vent, t_vent, v_vent

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_volcano_test_constructor(d)
  integer, intent(in) :: d

  integer :: fu, ierr
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.)   .or. &
       (mod_kinds_initialized.eqv..false.)      .or. &
       (mod_fu_manager_initialized.eqv..false.) .or. &
       (mod_physical_constants_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_volcano_test_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   ! Read input file
   call new_file_unit(fu,ierr)
   open(fu,file=trim(test_input_file_name), &
      status='old',action='read',form='formatted',iostat=ierr)
    if(ierr.ne.0) call error(this_sub_name,this_mod_name, &
      'Problems opening the input file')
    read(fu,input)
   close(fu,iostat=ierr)

   write(test_description(1),'(a,i1)') &
     'Volcano test, dimension d = ',d
   write(test_description(2),'(a,e9.3,a)') &
     '  viscosity nu = ',nu,';'
   write(test_description(3),'(3(a,e9.3),a)') &
     '  vent data: p=',p_vent,', T=',t_vent,', v=',v_vent,'.'

   ! No stratification and modified gas constants
   phc%omega   = 0.0_wp
   phc%gravity = 0.0_wp
   phc%rgas  = 18.46_wp
   phc%gamma = 1.02_wp
   phc%cv    = phc%rgas/(phc%gamma-1.0_wp)
   phc%cp    = phc%gamma*phc%cv
   phc%kappa = phc%rgas/phc%cp

   t_ref = 298.0_wp

   mod_volcano_test_initialized = .true.
 end subroutine mod_volcano_test_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_volcano_test_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_volcano_test_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_volcano_test_initialized = .false.
 end subroutine mod_volcano_test_destructor

!-----------------------------------------------------------------------

  !> Viscosity
  pure function coeff_visc(x,uu) result(mu)
   real(wp), intent(in) :: x(:,:), uu(:,:)
   real(wp) :: mu(2,size(x,2))

    mu = nu
  end function coeff_visc

!-----------------------------------------------------------------------

 !> Uniform and rest atmosphere
 pure function coeff_init(x,u_r) result(u0)
  real(wp), intent(in) :: x(:,:), u_r(:,:)
  real(wp) :: u0(2+size(x,1)+ntrcs,size(x,2))

   u0(1,:) = phc%p_s/(phc%rgas*t_ref) ! rho
   u0(2,:) = u0(1,:)*phc%cv*t_ref     ! e
   u0(3:,:) = 0.0_wp                  ! Ux

   ! subtract the reference state
   u0(1,:) = u0(1,:) - u_r(1,:) ! rho
   u0(2,:) = u0(2,:) - u_r(2,:) ! e

 end function coeff_init
 
!-----------------------------------------------------------------------

 !> Use the initial condition.
 pure function coeff_outref(x,u_r) result(u0)
  real(wp), intent(in) :: x(:,:), u_r(:,:)
  real(wp) :: u0(2+size(x,1)+ntrcs,size(x,2))

   u0 = coeff_init(x,u_r)

 end function coeff_outref

!-----------------------------------------------------------------------

 !> Dirichelt conditions at the vent
 pure function coeff_dir(x,u_r,breg) result(ub)
  real(wp), intent(in) :: x(:,:)
  real(wp), intent(in), optional :: u_r(:,:)
  integer, intent(in), optional :: breg
  real(wp) :: ub(2+size(x,1)+ntrcs,size(x,2))

   if(breg.eq.2) then ! this is the vent

     ! Overpressure and imposed velocity
     ub(1,:) = p_vent/(phc%rgas*t_vent)                    ! rho
     ub(2,:) = ub(1,:)*( phc%cv*t_ref + 0.5_wp*v_vent**2 ) ! e
     ub(3 ,:) = 0.0_wp                                     ! Ux
     ub(4 ,:) = ub(1,:)*v_vent                             ! Uy
     ub(5:,:) = 0.0_wp                                     ! Uz

     ! subtract the reference state
     ub(1,:) = ub(1,:) - u_r(1,:) ! rho
     ub(2,:) = ub(2,:) - u_r(2,:) ! e
   
   else
     ub = huge(ub) ! error
   endif

 end function coeff_dir

!-----------------------------------------------------------------------

 pure function coeff_neu(x,u_r,breg) result(h)
  real(wp), intent(in) :: x(:,:)
  real(wp), intent(in), optional :: u_r(:,:)
  integer, intent(in), optional :: breg
  real(wp) :: h(size(x,1),1+size(x,1)+ntrcs,size(x,2))

   h = 0.0_wp

 end function coeff_neu

!-----------------------------------------------------------------------

 !pure subroutine coeff_sponge(tmask,tau,taup,x)
 ! real(wp), intent(in) :: x(:,:)
 ! logical, intent(out) :: tmask
 ! real(wp), intent(out) :: tau(:), taup(:)

 ! real(wp), parameter :: pi_trig = 3.141592653589793_wp, &
 !   big = huge(1.0_wp)
 ! integer :: l, i
 ! real(wp) :: z, zn, d_bnd, coeff

 !  ! initialization
 !  tau = 0.0_wp
 !  taup = 0.0_wp
 !  tmask = .false.

 !  do l=1,size(x,2)

 !    ! upper layer
 !    z = x(size(x,1),l)
 !    if(z.gt.zd) then ! upper layer nonzero

 !      tmask = .true.

 !      ! gravity wave damping
 !      zn = (z-zd)/(zt-zd) ! normalized vertical coordinate
 !      if(zn.le.0.5_wp) then
 !        tau(l) = tau(l) + alpha/2.0_wp*(1.0_wp-cos(zn*pi_trig))
 !      else
 !        tau(l) = tau(l) + alpha/2.0_wp*(1.0_wp+(zn-0.5_wp)*pi_trig)
 !      endif

 !      ! acoustic wave damping
 !      d_bnd = (zt-z)/dv ! normalization
 !      if(d_bnd.le.3.0_wp) then
 !        ! tmask already set to .true.
 !        if(d_bnd.lt.1000.0_wp/big) then
 !          coeff = big/1000.0_wp
 !        else
 !          coeff = (1.0_wp-tanh(d_bnd))/tanh(d_bnd)
 !        endif
 !        tau(l)  = tau(l)  + coeff
 !        taup(l) = taup(l) + coeff
 !      endif

 !    endif

 !    ! lateral sponges (only acoustic wave damping)
 !    d_bnd = big
 !    do i=1,size(x,1)-1 ! first d-1 coordinates
 !      d_bnd = min( d_bnd , min( abs(x(i,l)-x1(i)) , abs(x(i,l)-x2(i)) ) )
 !    enddo
 !    d_bnd = d_bnd/dh ! normalization
 !    if(d_bnd.le.3.0_wp) then
 !      tmask = .true.
 !      if(d_bnd.lt.1000.0_wp/big) then
 !        coeff = big/1000.0_wp
 !      else
 !        coeff = (1.0_wp-tanh(d_bnd))/tanh(d_bnd)
 !      endif
 !      tau(l)  = tau(l)  + coeff
 !      taup(l) = taup(l) + coeff
 !    endif

 !  enddo
 !  
 !end subroutine coeff_sponge

!-----------------------------------------------------------------------

end module mod_volcano_test

